/* Bessel K_1(x) function in double precision */

#include <math.h>

double dbesk1(double x)
{
    int k;
    double t, y, v;
    static double a[16] = {
        1.5151605362537935201e-13, 3.363790951353651035e-11, 
        5.6514041131016827202e-9, 6.7816840255069534052e-7, 
        5.4253472222259226487e-5, 0.0026041666666666637057, 
        0.06250000000000000009, 0.5, 
        -8.9790303384748696588e-11, -1.4029047449249185771e-8, 
        -1.5592893752540998113e-6, -1.1253607018469017569e-4, 
        -0.0046421827665011579173, -0.085370719728648409609, 
        -0.3079657578292062966, 1.0000000000000000004
    };
    static double b[120] = {
        -9.4055461896630579928e-12, 3.1307934665844664773e-11, 
        4.2005295001519243251e-10, -4.1636196779679820012e-9, 
        1.4483026181700966164e-8, 1.1661000205428816914e-8, 
        -3.5023996724943046209e-7, 1.4404279316339005012e-6, 
        5.358156415715824208e-7, -3.5249754038612334639e-5, 
        1.7150324075631641453e-4, -4.1276362611239191024e-5, 
        -0.0046943110979636602591, 0.035085369853392357659, 
        0.20063574339907819159, 
        3.3998989888944034586e-11, 7.1558979072937373055e-11, 
        -2.9226856932927698732e-9, 1.4591620256525610213e-8, 
        -6.6141635609854161666e-9, -1.9991101838984472332e-7, 
        5.9185836628873530572e-7, 1.9562880347358085687e-6, 
        -1.5814366450418102764e-5, 7.6791682910944612028e-6, 
        2.8354678948323983936e-4, -0.0010217932669418690641, 
        -0.0032205661865210048433, 0.043497494842354644077, 
        0.16110284302315089935, 
        -2.0933987679665541827e-10, 7.9503322090520447316e-10, 
        3.8000150948242575774e-9, -2.3076136195585571309e-8, 
        -2.3902880302550799653e-8, 3.1914500937804377478e-7, 
        3.2639909831082417694e-7, -5.3166994792995439449e-6, 
        -3.1109524694269240094e-6, 9.2575906966353273247e-5, 
        7.5063709094147644746e-7, -0.0017416491592625765379, 
        0.0012138560335171676007, 0.045879687144659643175, 
        0.11566544716132846709, 
        3.1582384905164908749e-10, -1.9959561818098999516e-9, 
        8.6959328920030927557e-10, 1.1642778282445577109e-8, 
        4.3552264337818440471e-8, -1.5057982160481803238e-7, 
        -1.0101729117980989857e-6, 7.7002002510177612013e-7, 
        1.9580574235590194233e-5, 1.9358461980242834361e-5, 
        -3.3932339942485532728e-4, -9.3416673584325090073e-4, 
        0.0055800080455912847227, 0.038668683062477179235, 
        0.072651643500517000658, 
        -1.1554749629758510059e-10, 8.2270678758893273006e-10, 
        -5.0211156951551538591e-10, -1.4929179050554858361e-9, 
        -2.7107940791526366702e-8, -4.2204764086705349384e-8, 
        3.7253098167927628867e-7, 2.4374697215363361156e-6, 
        1.414194200690976837e-6, -4.8766389019473918231e-5, 
        -2.1681387247526720978e-4, 2.9325729929653405236e-4, 
        0.0064087534504827239815, 0.026054628289709454356, 
        0.040156431128194184336, 
        2.5506555170746221691e-11, -1.3521164018407978152e-10, 
        -8.3281235274106699399e-11, -9.7764849575562351891e-10, 
        3.4661828409940354542e-9, 3.9760633711791357544e-8, 
        1.590290664550452993e-7, -1.4919441249454941275e-7, 
        -5.3779684992094351263e-6, -2.7513862296246223142e-5, 
        -9.7880089725297162007e-6, 7.0787668964515789714e-4, 
        0.0046968199862345387583, 0.014745740181663320127, 
        0.020048622219583455723, 
        -3.4824483072529265585e-12, 1.5157161810563380451e-12, 
        8.5303859696700686144e-12, 3.3455414203743741076e-10, 
        2.0226016353844285376e-9, 5.312815400326633499e-9, 
        -3.0799322316418042137e-8, -4.4455408272954712128e-7, 
        -2.4293274626893384034e-6, -3.2129079340119038354e-6, 
        5.922540368307538885e-5, 5.6822962576781683532e-4, 
        0.0027152446516406682732, 0.0074075873691848838485, 
        0.0093044450815739269849, 
        -2.7683216166276377232e-13, 3.1986676777610155465e-12, 
        9.4142986954031445666e-12, 6.7934609179456399334e-11, 
        3.485152941147002933e-11, -2.5785248508896551557e-9, 
        -2.8310220027112571258e-8, -1.6384131113072271115e-7, 
        -3.2521663350596379097e-7, 4.038138875762230716e-6, 
        5.1917606978077281001e-5, 3.3420027947470126154e-4, 
        0.0013699550623118247094, 0.0034405619148342271096, 
        0.0041042919106665762794
    };
    static double c[120] = {
        4.5281968025889407937e-12, 1.0806749918195271176e-11, 
        9.6200972728717669027e-11, 5.721422706362526365e-10, 
        3.6077804282954825099e-9, 2.2465236858536681852e-8, 
        1.367696126430873523e-7, 7.9561767489531997361e-7, 
        4.3014380065615550573e-6, 2.092171390555028559e-5, 
        8.8079183950590176926e-5, 3.0549414408830252064e-4, 
        8.1295715613927890473e-4, 0.0014679809476357079195, 
        0.0013439197177355090057, 
        7.6019964430402432637e-13, -2.261619859915827119e-13, 
        1.7904450823779000744e-11, 9.1467054855312232717e-11, 
        7.1378582044879519122e-10, 4.9925255415445769102e-9, 
        3.3767315471315546644e-8, 2.1350774539167751457e-7, 
        1.2314353082655232903e-6, 6.2918685053670619181e-6, 
        2.7493229298777000013e-5, 9.8085825401369821771e-5, 
        2.6670282677770444935e-4, 4.8967895428135985381e-4, 
        4.5418248688489697144e-4, 
        9.4180115230375147213e-14, 7.5943117003734061145e-14, 
        3.0335730243874287654e-12, 2.0202796115462268051e-11, 
        1.6839020189186971198e-10, 1.2907875663127201526e-9, 
        9.354767612586579892e-9, 6.2471974110281880722e-8, 
        3.7585985422997380441e-7, 1.9838348288114906484e-6, 
        8.8884862203671982034e-6, 3.2333259238682810218e-5, 
        8.9266668913380400243e-5, 1.6589185669844051903e-4, 
        1.5536921180500113394e-4, 
        1.5425475332301107271e-14, 2.8674534590132490434e-14, 
        6.5078462279160216936e-13, 5.0939757793961391211e-12, 
        4.497983746074897552e-11, 3.6662925847520171711e-10, 
        2.7848878755089582413e-9, 1.929812005933947782e-8, 
        1.1950323861976892013e-7, 6.4513432758147478287e-7, 
        2.9422095033982461936e-6, 1.0854433321174584937e-5, 
        3.0307433185818899481e-5, 5.684098144306501785e-5, 
        5.3637016379451945253e-5, 
        3.1077953698439839352e-15, 8.6899496170729520378e-15, 
        1.6258562067326054104e-13, 1.4104842571366761537e-12, 
        1.3019455544084110747e-11, 1.1070466372863950239e-10, 
        8.6890603844230597917e-10, 6.1793722175049967488e-9, 
        3.9058865943755615801e-8, 2.1432806981070368523e-7, 
        9.9034657762983230155e-7, 3.6925185861895664251e-6, 
        1.0399877577259449786e-5, 1.9644939661550210015e-5, 
        1.8648773453825584597e-5, 
        7.2831555285336286457e-16, 2.6077534095895783532e-15, 
        4.4881202059263153495e-14, 4.1674329383944385626e-13, 
        3.9760100480223728037e-12, 3.483597635535118301e-11, 
        2.79932542127702497e-10, 2.0286513276830758107e-9, 
        1.3018343087118439152e-8, 7.2315927974997999365e-8, 
        3.3750708681924201599e-7, 1.2688020879407355571e-6, 
        3.5980954090811587848e-6, 6.8358260635246667316e-6, 
        6.5208606745808860557e-6, 
        1.9026412343503745875e-16, 8.0073765508732553766e-16, 
        1.3245754278055523992e-14, 1.2885201653055058502e-13, 
        1.2600129301230402587e-12, 1.1283306843147549277e-11, 
        9.2261481309646814329e-11, 6.7812033168299846818e-10, 
        4.4020645304595102132e-9, 2.4685719238301517679e-8, 
        1.1611886719473112509e-7, 4.3940380936523135466e-7, 
        1.2529878285546791905e-6, 2.3917218527087570384e-6, 
        2.290757464767187816e-6, 
        5.3709522135744366512e-17, 2.5239221050372845433e-16, 
        4.093314714589908336e-15, 4.1152784247617592367e-14, 
        4.0998840572769381012e-13, 3.7319354625807158852e-12, 
        3.0921671702920868014e-11, 2.2975898538634445343e-10, 
        1.5049754445782364328e-9, 8.5030864719789148982e-9, 
        4.025055939111842381e-8, 1.5312755642491878591e-7, 
        4.3865020375297892208e-7, 8.4059737392822153101e-7, 
        8.0785884122023473319e-7
    };
    static double d[40] = {
        9.2371554649979581914e-14, -2.3111336195788410887e-12, 
        5.7728710326649832559e-11, -1.8002298130091372598e-9, 
        7.6810375010517145638e-8, -5.2669973752193823306e-6, 
        0.0010112634961227401357, 0.16180215937964160466, 
        6.1381146507252683381e-14, -2.1034499679806301862e-12, 
        5.7090233460448415278e-11, -1.7990724350642330817e-9, 
        7.6809056078388019946e-8, -5.2669964425290062357e-6, 
        0.001011263495747828339, 0.16180215937970716383, 
        4.2458150578401296419e-14, -1.8435733128339016981e-12, 
        5.5534955081564656595e-11, -1.7938162188526358466e-9, 
        7.6798230945934117807e-8, -5.2669828728791921259e-6, 
        0.0010112634861753356559, 0.16180215938263409582, 
        3.0314798962267007518e-14, -1.5915009905364214455e-12, 
        5.3275907427402047438e-11, -1.7824862013841369751e-9, 
        7.676389035644707581e-8, -5.2669199860465945909e-6, 
        0.0010112634217687349189, 0.16180215941108227283, 
        2.2211515002229271212e-14, -1.3664088221521734796e-12, 
        5.0585177270502341602e-11, -1.7645432205894533462e-9, 
        7.6691805594577373698e-8, -5.2667455286976269634e-6, 
        0.001011263186281097458, 0.16180215954783127877
    };

    if (x < 0.8) {
        t = x * x;
        y = (((((((a[0] * t + a[1]) * t + 
            a[2]) * t + a[3]) * t + a[4]) * t + 
            a[5]) * t + a[6]) * t + a[7]) * x;
        y = (((((((a[8] * t + a[9]) * t + 
            a[10]) * t + a[11]) * t + a[12]) * t + 
            a[13]) * t + a[14]) * t + a[15]) / x + 
            y * log(x);
    } else if (x < 5.5) {
        v = 3 / x;
        t = x - v;
        k = (int) (t + 3);
        t = (k - 2) - t;
        k *= 15;
        y = ((((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * t + b[k + 13]) * t + 
            b[k + 14]) * v;
    } else if (x < 12.5) {
        k = (int) x;
        t = (k + 1) - x;
        k = 15 * (k - 5);
        y = (((((((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + 
            c[k + 11]) * t + c[k + 12]) * t + c[k + 13]) * t + 
            c[k + 14];
    } else {
        t = 60 / x;
        k = 8 * ((int) t);
        y = (((((((d[k] * t + d[k + 1]) * t + 
            d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + 
            d[k + 5]) * t + d[k + 6]) * t + d[k + 7]) * 
            sqrt(t) * exp(-x);
    }
    return y;
}

